#include <cmath>
#include <vector>

#include "macros.hpp"
#include "mdvector.hpp"
#include "points.hpp"

// TODO: Can get open-source scripts for generating quadrature points for arbitrary P>0
std::vector<double> Gauss_Legendre_pts(unsigned int P)
{
  std::vector<double> zeros(P,0.0);

  switch(P)
  {
    case 0:
      break; // do nothing
    case 1:
      zeros = {0.0}; break;

    case 2:
      zeros = {-std::sqrt(1./3.),
                std::sqrt(1./3.)}; break;

    case 3:
      zeros = {-std::sqrt(3./5.), 
               0.0, 
               std::sqrt(3./5.)}; break;

    case 4:
      zeros = {-std::sqrt((3.+2.*std::sqrt(6./5.))/7.), 
               -std::sqrt((3.-2.*std::sqrt(6./5.))/7.),
                std::sqrt((3.-2.*std::sqrt(6./5.))/7.), 
                std::sqrt((3.+2.*std::sqrt(6./5.))/7.)}; break;

    case 5:
      zeros = {-(1./3.)*std::sqrt(5.+2.*std::sqrt(10./7.)), 
               -(1./3.)*std::sqrt(5.-2.*std::sqrt(10./7.)),
                0.0,
                (1./3.)*std::sqrt(5.-2.*std::sqrt(10./7.)), 
                (1./3.)*std::sqrt(5.+2.*std::sqrt(10./7.))}; break;

    case 6:
      zeros = { -0.932469514203152,
                -0.661209386466264,
                -0.238619186083197,
                0.238619186083197,
                0.661209386466264,
                0.932469514203152}; break;
      
    case 7:
      zeros = { -0.949107912342758,
                -0.741531185599394,
                -0.405845151377397,
                0,
                0.405845151377397,
                0.741531185599394,
                0.949107912342758}; break;
      
    case 8:
      zeros = { -0.960289856497536,
                -0.796666477413627,
                -0.525532409916329,
                -0.183434642495650,
                0.183434642495650,
                0.525532409916329,
                0.796666477413627,
                0.960289856497536}; break;

    case 9:
      zeros = { -0.968160239507626,
                -0.836031107326636,
                -0.613371432700590,
                -0.324253423403809,
                0.0,
                0.324253423403809,
                0.613371432700590,
                0.836031107326636,    
                0.968160239507626}; break;

    case 10:
      zeros = {-0.973906528517172, 
               -0.865063366688985, 
               -0.679409568299024, 
               -0.433895394129247, 
               -0.148874338981631, 
               0.148874338981631,
               0.433895394129247, 
               0.679409568299024, 
               0.865063366688985, 
               0.973906528517172}; break;

    default:
      ThrowException("Gauss_Legendre_pts supports P up to 10!");
    
  }

  return zeros;
}

std::vector<double> DFRsp_pts(unsigned int P, double z1)
{
  std::vector<double> zeros(P,0.0);

  switch(P)
  {
    case 4:
      zeros = {-std::sqrt((3. - 5. * z1*z1)/(5. - 15.* z1 * z1)), 
               -z1,
                z1, 
                std::sqrt((3. - 5. * z1*z1)/(5. - 15.* z1 * z1))};
      break;

    default:
      ThrowException("DFRsp_pts supports limited cases right now!");
    
  }

  return zeros;

}

// TODO: Can get open-source scripts for generating quadrature weights for arbitrary P>0
std::vector<double> Gauss_Legendre_weights(unsigned int n)
{
  std::vector<double> weights(n,0.0);

  switch(n)
  {
    case 0:
      break; // do nothing
    case 1:
      weights = {2.0}; break;

    case 2:
      weights = {1.0, 1.0}; break;

    case 3:
      weights = {5./9., 8./9., 5./9.}; break;

    case 4:
      weights = {1./36. * (18. - std::sqrt(30.)),
                 1./36. * (18. + std::sqrt(30.)),
                 1./36. * (18. + std::sqrt(30.)),
                 1./36. * (18. - std::sqrt(30.))}; break;

    case 5:
      weights = {1./900. * (322 - 13. * std::sqrt(70.)),
                 1./900. * (322 + 13. * std::sqrt(70.)),
                 128./225.,
                 1./900. * (322 + 13. * std::sqrt(70.)),
                 1./900. * (322 - 13. * std::sqrt(70.))}; break;

    case 6:
      weights = {0.171324492379171,
                 0.360761573048139,
                 0.467913934572691,
                 0.467913934572691,
                 0.360761573048139,
                 0.171324492379171}; break;
      
    case 7:
      weights = {0.129484966168870,
                 0.279705391489277,
                 0.381830050505119,
                 0.417959183673469,
                 0.381830050505119,
                 0.279705391489277,
                 0.129484966168870}; break;
      
    case 8:
      weights = {0.101228536290377,
                 0.222381034453374,
                 0.313706645877887,
                 0.362683783378362,
                 0.362683783378362,
                 0.313706645877887,
                 0.222381034453374,
                 0.101228536290377}; break;
      
    case 9:
      weights = {0.081274388361575,
                 0.180648160694857,
                 0.260610696402936,
                 0.312347077040003,
                 0.330239355001260,
                 0.312347077040003,
                 0.260610696402936,
                 0.180648160694857,
                 0.081274388361575}; break;
      
    case 10:
      weights = {0.066671344308688, 
                 0.149451349150581,
                 0.219086362515982,
                 0.269266719309996, 
                 0.295524224714753, 
                 0.295524224714753, 
                 0.269266719309996, 
                 0.219086362515982, 
                 0.149451349150581,
                 0.066671344308688}; break;

    default:
      ThrowException("Gauss_Legendre_weights supports up to 10 points!");

  }

  return weights;
}

std::vector<double> Shape_pts(unsigned int P)
{
  std::vector<double> nodes(P+1, 0.0);

  double dx = 2.0/(double)P;

  nodes[0] = -1.0;
  for (unsigned int i = 1; i < P; i++)
    nodes[i] = nodes[i-1] + dx;

  nodes[P] = 1.0;

  return nodes; 
}

mdvector<double> WS_Tri_pts(unsigned int P)
{
  unsigned int nPts = (P+1)*(P+2)/2;
  mdvector<double> pts({nPts, 2});

  switch(P)
  {
    case 0:
      pts(0, 0) = -0.333333333333; pts(0, 1) = -0.333333333333;
      break;

    case 1:
      pts( 0, 0) = -0.6666666666666660; pts( 0, 1) =  0.3333333333333320; 
      pts( 1, 0) =  0.3333333333333320; pts( 1, 1) = -0.6666666666666660; 
      pts( 2, 0) = -0.6666666666666660; pts( 2, 1) = -0.6666666666666660; 
      break;

    case 2:
      pts( 0, 0) = -0.8168475729804400; pts( 0, 1) =  0.6336951459608799; 
      pts( 1, 0) =  0.6336951459608799; pts( 1, 1) = -0.8168475729804400; 
      pts( 2, 0) = -0.8168475729804400; pts( 2, 1) = -0.8168475729804400; 
      pts( 3, 0) = -0.1081030181680720; pts( 3, 1) = -0.7837939636638560; 
      pts( 4, 0) = -0.7837939636638560; pts( 4, 1) = -0.1081030181680720; 
      pts( 5, 0) = -0.1081030181680720; pts( 5, 1) = -0.1081030181680720; 
      break;

    case 3:
      pts( 0, 0) = -0.3333333333333330; pts( 0, 1) = -0.3333333333333330; 
      pts( 1, 0) = -0.8888718946604141; pts( 1, 1) =  0.7777437893208280; 
      pts( 2, 0) =  0.7777437893208280; pts( 2, 1) = -0.8888718946604141; 
      pts( 3, 0) = -0.8888718946604141; pts( 3, 1) = -0.8888718946604141; 
      pts( 4, 0) =  0.2684214954914460; pts( 4, 1) = -0.4089325765282140; 
      pts( 5, 0) = -0.8594889189632320; pts( 5, 1) = -0.4089325765282140; 
      pts( 6, 0) = -0.4089325765282140; pts( 6, 1) = -0.8594889189632320; 
      pts( 7, 0) =  0.2684214954914460; pts( 7, 1) = -0.8594889189632320; 
      pts( 8, 0) = -0.8594889189632320; pts( 8, 1) =  0.2684214954914460; 
      pts( 9, 0) = -0.4089325765282140; pts( 9, 1) =  0.2684214954914460; 
      break;

    case 4:
      pts( 0, 0) = -0.9282582446085320; pts( 0, 1) =  0.8565164892170640; 
      pts( 1, 0) =  0.8565164892170640; pts( 1, 1) = -0.9282582446085320; 
      pts( 2, 0) = -0.9282582446085320; pts( 2, 1) = -0.9282582446085320; 
      pts( 3, 0) = -0.5165412084640660; pts( 3, 1) =  0.0330824169281320; 
      pts( 4, 0) =  0.0330824169281320; pts( 4, 1) = -0.5165412084640660; 
      pts( 5, 0) = -0.5165412084640660; pts( 5, 1) = -0.5165412084640660; 
      pts( 6, 0) = -0.0513824244458420; pts( 6, 1) = -0.8972351511083160; 
      pts( 7, 0) = -0.8972351511083160; pts( 7, 1) = -0.0513824244458420; 
      pts( 8, 0) = -0.0513824244458420; pts( 8, 1) = -0.0513824244458420; 
      pts( 9, 0) =  0.5023672622129680; pts( 9, 1) = -0.5969922362364000; 
      pts(10, 0) = -0.9053750259765680; pts(10, 1) = -0.5969922362364000; 
      pts(11, 0) = -0.5969922362364000; pts(11, 1) = -0.9053750259765680; 
      pts(12, 0) =  0.5023672622129680; pts(12, 1) = -0.9053750259765680; 
      pts(13, 0) = -0.9053750259765680; pts(13, 1) =  0.5023672622129680; 
      pts(14, 0) = -0.5969922362364000; pts(14, 1) =  0.5023672622129680; 
      break;

    case 5:
      pts( 0, 0) = -0.9437740956346720; pts( 0, 1) =  0.8875481912693440; 
      pts( 1, 0) =  0.8875481912693440; pts( 1, 1) = -0.9437740956346720; 
      pts( 2, 0) = -0.9437740956346720; pts( 2, 1) = -0.9437740956346720; 
      pts( 3, 0) = -0.6457218030613660; pts( 3, 1) =  0.2914436061227320; 
      pts( 4, 0) =  0.2914436061227320; pts( 4, 1) = -0.6457218030613660; 
      pts( 5, 0) = -0.6457218030613660; pts( 5, 1) = -0.6457218030613660; 
      pts( 6, 0) = -0.1889828082651340; pts( 6, 1) = -0.6220343834697321; 
      pts( 7, 0) = -0.6220343834697321; pts( 7, 1) = -0.1889828082651340; 
      pts( 8, 0) = -0.1889828082651340; pts( 8, 1) = -0.1889828082651340; 
      pts( 9, 0) =  0.6358019600569980; pts( 9, 1) = -0.7028683754582260; 
      pts(10, 0) = -0.9329335845987720; pts(10, 1) = -0.7028683754582260; 
      pts(11, 0) = -0.7028683754582260; pts(11, 1) = -0.9329335845987720; 
      pts(12, 0) =  0.6358019600569980; pts(12, 1) = -0.9329335845987720; 
      pts(13, 0) = -0.9329335845987720; pts(13, 1) =  0.6358019600569980; 
      pts(14, 0) = -0.7028683754582260; pts(14, 1) =  0.6358019600569980; 
      pts(15, 0) =  0.2099578235502660; pts(15, 1) = -0.2856074027686380; 
      pts(16, 0) = -0.9243504207816280; pts(16, 1) = -0.2856074027686380; 
      pts(17, 0) = -0.2856074027686380; pts(17, 1) = -0.9243504207816280; 
      pts(18, 0) =  0.2099578235502660; pts(18, 1) = -0.9243504207816280; 
      pts(19, 0) = -0.9243504207816280; pts(19, 1) =  0.2099578235502660; 
      pts(20, 0) = -0.2856074027686380; pts(20, 1) =  0.2099578235502660; 
      break;

    case 6:
      pts( 0, 0) = -0.3333333333333330; pts( 0, 1) = -0.3333333333333330; 
      pts( 1, 0) = -0.9600456257556140; pts( 1, 1) =  0.9200912515112279; 
      pts( 2, 0) =  0.9200912515112279; pts( 2, 1) = -0.9600456257556140; 
      pts( 3, 0) = -0.9600456257556140; pts( 3, 1) = -0.9600456257556140; 
      pts( 4, 0) = -0.7365564649400040; pts( 4, 1) =  0.4731129298800080; 
      pts( 5, 0) =  0.4731129298800080; pts( 5, 1) = -0.7365564649400040; 
      pts( 6, 0) = -0.7365564649400040; pts( 6, 1) = -0.7365564649400040; 
      pts( 7, 0) = -0.0297293064130780; pts( 7, 1) = -0.9405413871738439; 
      pts( 8, 0) = -0.9405413871738439; pts( 8, 1) = -0.0297293064130780; 
      pts( 9, 0) = -0.0297293064130780; pts( 9, 1) = -0.0297293064130780; 
      pts(10, 0) =  0.7358224202359001; pts(10, 1) = -0.7840960363079780; 
      pts(11, 0) = -0.9517263839279220; pts(11, 1) = -0.7840960363079780; 
      pts(12, 0) = -0.7840960363079780; pts(12, 1) = -0.9517263839279220; 
      pts(13, 0) =  0.7358224202359001; pts(13, 1) = -0.9517263839279220; 
      pts(14, 0) = -0.9517263839279220; pts(14, 1) =  0.7358224202359001; 
      pts(15, 0) = -0.7840960363079780; pts(15, 1) =  0.7358224202359001; 
      pts(16, 0) =  0.4017451407614460; pts(16, 1) = -0.4583184541568660; 
      pts(17, 0) = -0.9434266866045800; pts(17, 1) = -0.4583184541568660; 
      pts(18, 0) = -0.4583184541568660; pts(18, 1) = -0.9434266866045800; 
      pts(19, 0) =  0.4017451407614460; pts(19, 1) = -0.9434266866045800; 
      pts(20, 0) = -0.9434266866045800; pts(20, 1) =  0.4017451407614460; 
      pts(21, 0) = -0.4583184541568660; pts(21, 1) =  0.4017451407614460; 
      pts(22, 0) =  0.0733093684122760; pts(22, 1) = -0.3669008023107660; 
      pts(23, 0) = -0.7064085661015100; pts(23, 1) = -0.3669008023107660; 
      pts(24, 0) = -0.3669008023107660; pts(24, 1) = -0.7064085661015100; 
      pts(25, 0) =  0.0733093684122760; pts(25, 1) = -0.7064085661015100; 
      pts(26, 0) = -0.7064085661015100; pts(26, 1) =  0.0733093684122760; 
      pts(27, 0) = -0.3669008023107660; pts(27, 1) =  0.0733093684122760; 
      break;

    case 7:
      pts( 0, 0) = -0.9576571544410700; pts( 0, 1) =  0.9153143088821400; 
      pts( 1, 0) =  0.9153143088821400; pts( 1, 1) = -0.9576571544410700; 
      pts( 2, 0) = -0.9576571544410700; pts( 2, 1) = -0.9576571544410700; 
      pts( 3, 0) = -0.7988312052082240; pts( 3, 1) =  0.5976624104164480; 
      pts( 4, 0) =  0.5976624104164480; pts( 4, 1) = -0.7988312052082240; 
      pts( 5, 0) = -0.7988312052082240; pts( 5, 1) = -0.7988312052082240; 
      pts( 6, 0) = -0.4579233845761360; pts( 6, 1) = -0.0841532308477280; 
      pts( 7, 0) = -0.0841532308477280; pts( 7, 1) = -0.4579233845761360; 
      pts( 8, 0) = -0.4579233845761360; pts( 8, 1) = -0.4579233845761360; 
      pts( 9, 0) = -0.1196174831923360; pts( 9, 1) = -0.7607650336153280; 
      pts(10, 0) = -0.7607650336153280; pts(10, 1) = -0.1196174831923360; 
      pts(11, 0) = -0.1196174831923360; pts(11, 1) = -0.1196174831923360; 
      pts(12, 0) =  0.7599592828544620; pts(12, 1) = -0.9634866418505040; 
      pts(13, 0) = -0.7964726410039580; pts(13, 1) = -0.9634866418505040; 
      pts(14, 0) = -0.9634866418505040; pts(14, 1) = -0.7964726410039580; 
      pts(15, 0) =  0.7599592828544620; pts(15, 1) = -0.7964726410039580; 
      pts(16, 0) = -0.7964726410039580; pts(16, 1) =  0.7599592828544620; 
      pts(17, 0) = -0.9634866418505040; pts(17, 1) =  0.7599592828544620; 
      pts(18, 0) =  0.1651240457273440; pts(18, 1) = -0.2119334566600260; 
      pts(19, 0) = -0.9531905890673180; pts(19, 1) = -0.2119334566600260; 
      pts(20, 0) = -0.2119334566600260; pts(20, 1) = -0.9531905890673180; 
      pts(21, 0) =  0.1651240457273440; pts(21, 1) = -0.9531905890673180; 
      pts(22, 0) = -0.9531905890673180; pts(22, 1) =  0.1651240457273440; 
      pts(23, 0) = -0.2119334566600260; pts(23, 1) =  0.1651240457273440; 
      pts(24, 0) =  0.5030612290855641; pts(24, 1) = -0.5475089381815420; 
      pts(25, 0) = -0.9555522909040221; pts(25, 1) = -0.5475089381815420; 
      pts(26, 0) = -0.5475089381815420; pts(26, 1) = -0.9555522909040221; 
      pts(27, 0) =  0.5030612290855641; pts(27, 1) = -0.9555522909040221; 
      pts(28, 0) = -0.9555522909040221; pts(28, 1) =  0.5030612290855641; 
      pts(29, 0) = -0.5475089381815420; pts(29, 1) =  0.5030612290855641; 
      pts(30, 0) =  0.2714743665262100; pts(30, 1) = -0.5018415447573360; 
      pts(31, 0) = -0.7696328217688740; pts(31, 1) = -0.5018415447573360; 
      pts(32, 0) = -0.5018415447573360; pts(32, 1) = -0.7696328217688740; 
      pts(33, 0) =  0.2714743665262100; pts(33, 1) = -0.7696328217688740; 
      pts(34, 0) = -0.7696328217688740; pts(34, 1) =  0.2714743665262100; 
      pts(35, 0) = -0.5018415447573360; pts(35, 1) =  0.2714743665262100; 
      break;

    case 8:
      pts( 0, 0) = -0.8958458673313620; pts( 0, 1) =  0.7916917346627240; 
      pts( 1, 0) =  0.7916917346627240; pts( 1, 1) = -0.8958458673313620; 
      pts( 2, 0) = -0.8958458673313620; pts( 2, 1) = -0.8958458673313620; 
      pts( 3, 0) = -0.7725869427780080; pts( 3, 1) =  0.5451738855560160; 
      pts( 4, 0) =  0.5451738855560160; pts( 4, 1) = -0.7725869427780080; 
      pts( 5, 0) = -0.7725869427780080; pts( 5, 1) = -0.7725869427780080; 
      pts( 6, 0) = -0.5656623084521640; pts( 6, 1) =  0.1313246169043280; 
      pts( 7, 0) =  0.1313246169043280; pts( 7, 1) = -0.5656623084521640; 
      pts( 8, 0) = -0.5656623084521640; pts( 8, 1) = -0.5656623084521640; 
      pts( 9, 0) = -0.2278544512773540; pts( 9, 1) = -0.5442910974452920; 
      pts(10, 0) = -0.5442910974452920; pts(10, 1) = -0.2278544512773540; 
      pts(11, 0) = -0.2278544512773540; pts(11, 1) = -0.2278544512773540; 
      pts(12, 0) = -0.0201481341734600; pts(12, 1) = -0.9597037316530800; 
      pts(13, 0) = -0.9597037316530800; pts(13, 1) = -0.0201481341734600; 
      pts(14, 0) = -0.0201481341734600; pts(14, 1) = -0.0201481341734600; 
      pts(15, 0) =  0.6693761812315560; pts(15, 1) = -0.7035037567771640; 
      pts(16, 0) = -0.9658724244543920; pts(16, 1) = -0.7035037567771640; 
      pts(17, 0) = -0.7035037567771640; pts(17, 1) = -0.9658724244543920; 
      pts(18, 0) =  0.6693761812315560; pts(18, 1) = -0.9658724244543920; 
      pts(19, 0) = -0.9658724244543920; pts(19, 1) =  0.6693761812315560; 
      pts(20, 0) = -0.7035037567771640; pts(20, 1) =  0.6693761812315560; 
      pts(21, 0) =  0.4144959853084960; pts(21, 1) = -0.5795767423228601; 
      pts(22, 0) = -0.8349192429856360; pts(22, 1) = -0.5795767423228601; 
      pts(23, 0) = -0.5795767423228601; pts(23, 1) = -0.8349192429856360; 
      pts(24, 0) =  0.4144959853084960; pts(24, 1) = -0.8349192429856360; 
      pts(25, 0) = -0.8349192429856360; pts(25, 1) =  0.4144959853084960; 
      pts(26, 0) = -0.5795767423228601; pts(26, 1) =  0.4144959853084960; 
      pts(27, 0) =  0.9258038332858580; pts(27, 1) = -0.9258815152162200; 
      pts(28, 0) = -0.9999223180696380; pts(28, 1) = -0.9258815152162200; 
      pts(29, 0) = -0.9258815152162200; pts(29, 1) = -0.9999223180696380; 
      pts(30, 0) =  0.9258038332858580; pts(30, 1) = -0.9999223180696380; 
      pts(31, 0) = -0.9999223180696380; pts(31, 1) =  0.9258038332858580; 
      pts(32, 0) = -0.9258815152162200; pts(32, 1) =  0.9258038332858580; 
      pts(33, 0) =  0.3480122228257960; pts(33, 1) = -0.3820797590136400; 
      pts(34, 0) = -0.9659324638121560; pts(34, 1) = -0.3820797590136400; 
      pts(35, 0) = -0.3820797590136400; pts(35, 1) = -0.9659324638121560; 
      pts(36, 0) =  0.3480122228257960; pts(36, 1) = -0.9659324638121560; 
      pts(37, 0) = -0.9659324638121560; pts(37, 1) =  0.3480122228257960; 
      pts(38, 0) = -0.3820797590136400; pts(38, 1) =  0.3480122228257960; 
      pts(39, 0) =  0.0838979176456600; pts(39, 1) = -0.2793516463213580; 
      pts(40, 0) = -0.8045462713243020; pts(40, 1) = -0.2793516463213580; 
      pts(41, 0) = -0.2793516463213580; pts(41, 1) = -0.8045462713243020; 
      pts(42, 0) =  0.0838979176456600; pts(42, 1) = -0.8045462713243020; 
      pts(43, 0) = -0.8045462713243020; pts(43, 1) =  0.0838979176456600; 
      pts(44, 0) = -0.2793516463213580; pts(44, 1) =  0.0838979176456600; 
      break;

    default:
      ThrowException("WS_Tri_pts undefined for provided order!");
  }

  return pts;
}

mdvector<double> WS_Tri_weights(unsigned int P)
{
  unsigned int nPts = (P+1)*(P+2)/2;
  mdvector<double> weights({nPts});

  switch(P)
  {
    case 0:
      weights(0) = 2.0;
      break;

    case 1:
      weights( 0) =  0.6666666666666660;
      weights( 1) =  0.6666666666666660;
      weights( 2) =  0.6666666666666660;
      break;

    case 2:
      weights( 0) =  0.2199034873106660;
      weights( 1) =  0.2199034873106660;
      weights( 2) =  0.2199034873106660;
      weights( 3) =  0.4467631793560000;
      weights( 4) =  0.4467631793560000;
      weights( 5) =  0.4467631793560000;
      break;

    case 3:
      weights( 0) =  0.4030859771694600;
      weights( 1) =  0.0839110259932980;
      weights( 2) =  0.0839110259932980;
      weights( 3) =  0.0839110259932980;
      weights( 4) =  0.2241968241417740;
      weights( 5) =  0.2241968241417740;
      weights( 6) =  0.2241968241417740;
      weights( 7) =  0.2241968241417740;
      weights( 8) =  0.2241968241417740;
      weights( 9) =  0.2241968241417740;
      break;

    case 4:
      weights( 0) =  0.0358309100246060;
      weights( 1) =  0.0358309100246060;
      weights( 2) =  0.0358309100246060;
      weights( 3) =  0.2554243917625300;
      weights( 4) =  0.2554243917625300;
      weights( 5) =  0.2554243917625300;
      weights( 6) =  0.1524121247710700;
      weights( 7) =  0.1524121247710700;
      weights( 8) =  0.1524121247710700;
      weights( 9) =  0.1114996200542300;
      weights(10) =  0.1114996200542300;
      weights(11) =  0.1114996200542300;
      weights(12) =  0.1114996200542300;
      weights(13) =  0.1114996200542300;
      weights(14) =  0.1114996200542300;
      break;

    case 5:
      weights( 0) =  0.0207187493930760;
      weights( 1) =  0.0207187493930760;
      weights( 2) =  0.0207187493930760;
      weights( 3) =  0.1507897686534760;
      weights( 4) =  0.1507897686534760;
      weights( 5) =  0.1507897686534760;
      weights( 6) =  0.1950956047464840;
      weights( 7) =  0.1950956047464840;
      weights( 8) =  0.1950956047464840;
      weights( 9) =  0.0579385387449460;
      weights(10) =  0.0579385387449460;
      weights(11) =  0.0579385387449460;
      weights(12) =  0.0579385387449460;
      weights(13) =  0.0579385387449460;
      weights(14) =  0.0579385387449460;
      weights(15) =  0.0920927331918700;
      weights(16) =  0.0920927331918700;
      weights(17) =  0.0920927331918700;
      weights(18) =  0.0920927331918700;
      weights(19) =  0.0920927331918700;
      weights(20) =  0.0920927331918700;
      break;

    case 6:
      weights( 0) =  0.1672164244312740;
      weights( 1) =  0.0105443405609900;
      weights( 2) =  0.0105443405609900;
      weights( 3) =  0.0105443405609900;
      weights( 4) =  0.0891058733590080;
      weights( 5) =  0.0891058733590080;
      weights( 6) =  0.0891058733590080;
      weights( 7) =  0.0676314256083960;
      weights( 8) =  0.0676314256083960;
      weights( 9) =  0.0676314256083960;
      weights(10) =  0.0314209226803660;
      weights(11) =  0.0314209226803660;
      weights(12) =  0.0314209226803660;
      weights(13) =  0.0314209226803660;
      weights(14) =  0.0314209226803660;
      weights(15) =  0.0314209226803660;
      weights(16) =  0.0564102725612320;
      weights(17) =  0.0564102725612320;
      weights(18) =  0.0564102725612320;
      weights(19) =  0.0564102725612320;
      weights(20) =  0.0564102725612320;
      weights(21) =  0.0564102725612320;
      weights(22) =  0.1339919142556600;
      weights(23) =  0.1339919142556600;
      weights(24) =  0.1339919142556600;
      weights(25) =  0.1339919142556600;
      weights(26) =  0.1339919142556600;
      weights(27) =  0.1339919142556600;
      break;

    case 7:
      weights( 0) =  0.0112782475738200;
      weights( 1) =  0.0112782475738200;
      weights( 2) =  0.0112782475738200;
      weights( 3) =  0.0542979363845560;
      weights( 4) =  0.0542979363845560;
      weights( 5) =  0.0542979363845560;
      weights( 6) =  0.1262018250667180;
      weights( 7) =  0.1262018250667180;
      weights( 8) =  0.1262018250667180;
      weights( 9) =  0.1035055913597980;
      weights(10) =  0.1035055913597980;
      weights(11) =  0.1035055913597980;
      weights(12) =  0.0197335071492920;
      weights(13) =  0.0197335071492920;
      weights(14) =  0.0197335071492920;
      weights(15) =  0.0197335071492920;
      weights(16) =  0.0197335071492920;
      weights(17) =  0.0197335071492920;
      weights(18) =  0.0440164096002940;
      weights(19) =  0.0440164096002940;
      weights(20) =  0.0440164096002940;
      weights(21) =  0.0440164096002940;
      weights(22) =  0.0440164096002940;
      weights(23) =  0.0440164096002940;
      weights(24) =  0.0332891401534720;
      weights(25) =  0.0332891401534720;
      weights(26) =  0.0332891401534720;
      weights(27) =  0.0332891401534720;
      weights(28) =  0.0332891401534720;
      weights(29) =  0.0332891401534720;
      weights(30) =  0.0886524762378280;
      weights(31) =  0.0886524762378280;
      weights(32) =  0.0886524762378280;
      weights(33) =  0.0886524762378280;
      weights(34) =  0.0886524762378280;
      weights(35) =  0.0886524762378280;
      break;

    case 8:
      weights( 0) =  0.0264519521410120;
      weights( 1) =  0.0264519521410120;
      weights( 2) =  0.0264519521410120;
      weights( 3) =  0.0347516292773620;
      weights( 4) =  0.0347516292773620;
      weights( 5) =  0.0347516292773620;
      weights( 6) =  0.0945461565676520;
      weights( 7) =  0.0945461565676520;
      weights( 8) =  0.0945461565676520;
      weights( 9) =  0.1060044409420360;
      weights(10) =  0.1060044409420360;
      weights(11) =  0.1060044409420360;
      weights(12) =  0.0378851639278180;
      weights(13) =  0.0378851639278180;
      weights(14) =  0.0378851639278180;
      weights(15) =  0.0238755661869220;
      weights(16) =  0.0238755661869220;
      weights(17) =  0.0238755661869220;
      weights(18) =  0.0238755661869220;
      weights(19) =  0.0238755661869220;
      weights(20) =  0.0238755661869220;
      weights(21) =  0.0493958430852700;
      weights(22) =  0.0493958430852700;
      weights(23) =  0.0493958430852700;
      weights(24) =  0.0493958430852700;
      weights(25) =  0.0493958430852700;
      weights(26) =  0.0493958430852700;
      weights(27) =  0.0043857672678220;
      weights(28) =  0.0043857672678220;
      weights(29) =  0.0043857672678220;
      weights(30) =  0.0043857672678220;
      weights(31) =  0.0043857672678220;
      weights(32) =  0.0043857672678220;
      weights(33) =  0.0308750382999060;
      weights(34) =  0.0308750382999060;
      weights(35) =  0.0308750382999060;
      weights(36) =  0.0308750382999060;
      weights(37) =  0.0308750382999060;
      weights(38) =  0.0308750382999060;
      weights(39) =  0.0749814470654740;
      weights(40) =  0.0749814470654740;
      weights(41) =  0.0749814470654740;
      weights(42) =  0.0749814470654740;
      weights(43) =  0.0749814470654740;
      weights(44) =  0.0749814470654740;
      break;

    default:
      ThrowException("WS_Tri_weights undefined for provided order!");
  }

  return weights;
}

mdvector<double> RW_Tri_pts(unsigned int P)
{
  unsigned int nPts = (P+1)*(P+2)/2;
  mdvector<double> pts({nPts, 2});

  switch(P)
  {
    case 1:
      pts(0, 0) = -0.666666666667;
      pts(1, 0) = 0.333333333333;
      pts(2, 0) = -0.666666666667;

      pts(0, 1) = 0.333333333333;
      pts(1, 1) = -0.666666666667;
      pts(2, 1) = -0.666666666667;
      break;

    case 2:
      pts(0, 0) = -0.81684757298;
      pts(1, 0) = 0.633695145961;
      pts(2, 0) = -0.81684757298;
      pts(3, 0) = -0.108103018168;
      pts(4, 0) = -0.783793963664;
      pts(5, 0) = -0.108103018168;

      pts(0, 1) = 0.633695145961;
      pts(1, 1) = -0.81684757298;
      pts(2, 1) = -0.81684757298;
      pts(3, 1) = -0.783793963664;
      pts(4, 1) = -0.108103018168;
      pts(5, 1) = -0.108103018168;
      break;

    case 3:
      // WS
      pts = symm_pts_tri({1./3.}, {0.055564052669793}, {0.295533711735893}, {0.070255540518384});

      // FW err opt
      //pts = symm_pts_tri({1./3.}, {.05275582964828789}, {.293034268156967841039435067323947406114147778}, {.06998504422734275906007349011531234});

      // fmincon eps 1e-8 tolcon 1e-9 ws start
      //pts = symm_pts_tri({1./3.}, {0.055744804500109}, {0.290384222060825}, {0.070841762329827});

      // fmincon eps 1e-10 tolcon 1e-11 ws start
      //pts = symm_pts_tri({1./3.}, {0.055758983558155}, {0.290285227512689}, {0.070857385399496});

      break;

    case 4:
      // WS
      pts = symm_pts_tri({}, {0.035870877695734,0.241729395767967,0.474308787777079}, {0.201503881881800}, {0.047312487011716});

      // fmincon eps 1e-8 tolcon 1e-9 ws start
      //pts = symm_pts_tri({}, {0.034997794563496, 0.242941826150418, 0.473422728963121}, {0.200163370171128}, {0.047078563918947});

      // fmincon eps 1e-10 tolcon 1e-11 ws start
      //pts = symm_pts_tri({}, {0.034681580220044, 0.243071555674492, 0.473372556704605}, {0.200039998995093}, {0.047293668511439});
      
      break;

    default:
      ThrowException("RW_Tri_pts undefined for provided order!");
  }

  return pts;
}

mdvector<double> symm_pts_tri(std::vector<double> a3, std::vector<double> a21, std::vector<double> a111, std::vector<double> b111)
{
  auto n3 = a3.size();
  auto n21 = a21.size();
  auto n111 = a111.size();

  unsigned int nPts = n3 + n21 * 3 + n111 * 6;

  mdvector<double> pts({nPts, 2}, 0);

  std::vector<std::vector<double>> trans(2);
  trans[0] = {-1, 1, -1};
  trans[1] = {-1, -1, 1};

  int pt = 0;

  // n3 orbits
  for (unsigned int i = 0; i < n3; i++)
  {
    std::vector<double> lam = {a3[i], a3[i], a3[i]}; 

    for (unsigned int j = 0; j < 2; j++)
    {
      for (unsigned int k = 0; k < 3; k++)
      {
        pts(pt, j) += lam[k] * trans[j][k];
      }
    }

    pt++;  
  }

  // n21 orbits
  std::vector<std::vector<int>> perm21(3);
  perm21[0] = {0, 1, 2};
  perm21[1] = {2, 0, 1};
  perm21[2] = {0, 2, 1};

  for (unsigned int i = 0; i < n21; i++)
  {
    std::vector<double> lam0 = {a21[i], a21[i], 1 - 2 * a21[i]}; 

    for (unsigned int n = 0; n < 3; n++)
    {
      auto lam = lam0;
      for (unsigned int m = 0; m < 3; m++)
        lam[m] = lam0[perm21[n][m]];

      for (unsigned int j = 0; j < 2; j++)
      {
        for (unsigned int k = 0; k < 3; k++)
        {
          pts(pt, j) += lam[k] * trans[j][k];
        }
      }

      pt++;  
    }
  }

  // n111 orbits
  std::vector<std::vector<int>> perm111(6);
  perm111[0] = {0, 1, 2};
  perm111[1] = {2, 0, 1};
  perm111[2] = {0, 2, 1};
  perm111[3] = {2, 1, 0};
  perm111[4] = {1, 0, 2};
  perm111[5] = {1, 2, 0};

  for (unsigned int i = 0; i < n111; i++)
  {
    std::vector<double> lam0 = {a111[i], b111[i], 1 - a111[i] - b111[i]}; 

    for (unsigned int n = 0; n < 6; n++)
    {
      auto lam = lam0;
      for (unsigned int m = 0; m < 3; m++)
        lam[m] = lam0[perm111[n][m]];

      for (unsigned int j = 0; j < 2; j++)
      {
        for (unsigned int k = 0; k < 3; k++)
        {
          pts(pt, j) += lam[k] * trans[j][k];
        }
      }

      pt++;  
    }
  }

  return pts;
  
}

mdvector<double> WS_Tet_pts(unsigned int P)
{
  unsigned int nPts = (P+1)*(P+2)*(P+3)/6;
  mdvector<double> pts({nPts, 3});

  switch(P)
  {
    case 0:
      pts(0, 0) = -0.5; pts(0, 1) = -0.5; pts(0, 2) = -0.5;
      break;

    case 1:
      pts( 0, 0) = -0.7236067977499780; pts( 0, 1) = -0.7236067977499780; pts( 0, 2) =  0.1708203932499340; 
      pts( 1, 0) =  0.1708203932499340; pts( 1, 1) = -0.7236067977499780; pts( 1, 2) = -0.7236067977499780; 
      pts( 2, 0) = -0.7236067977499780; pts( 2, 1) =  0.1708203932499340; pts( 2, 2) = -0.7236067977499780; 
      pts( 3, 0) = -0.7236067977499780; pts( 3, 1) = -0.7236067977499780; pts( 3, 2) = -0.7236067977499780; 
      break;

    case 2:
      pts( 0, 0) = -0.8523301965475532; pts( 0, 1) = -0.8523301965475532; pts( 0, 2) =  0.5569905896426596; 
      pts( 1, 0) =  0.5569905896426596; pts( 1, 1) = -0.8523301965475532; pts( 1, 2) = -0.8523301965475532; 
      pts( 2, 0) = -0.8523301965475532; pts( 2, 1) =  0.5569905896426596; pts( 2, 2) = -0.8523301965475532; 
      pts( 3, 0) = -0.8523301965475532; pts( 3, 1) = -0.8523301965475532; pts( 3, 2) = -0.8523301965475532; 
      pts( 4, 0) = -0.1875113122318980; pts( 4, 1) = -0.8124886877681020; pts( 4, 2) = -0.1875113122318980; 
      pts( 5, 0) = -0.8124886877681020; pts( 5, 1) = -0.1875113122318980; pts( 5, 2) = -0.1875113122318980; 
      pts( 6, 0) = -0.1875113122318980; pts( 6, 1) = -0.8124886877681020; pts( 6, 2) = -0.8124886877681020; 
      pts( 7, 0) = -0.8124886877681020; pts( 7, 1) = -0.1875113122318980; pts( 7, 2) = -0.8124886877681020; 
      pts( 8, 0) = -0.8124886877681020; pts( 8, 1) = -0.8124886877681020; pts( 8, 2) = -0.1875113122318980; 
      pts( 9, 0) = -0.1875113122318980; pts( 9, 1) = -0.1875113122318980; pts( 9, 2) = -0.8124886877681020; 
      break;

    case 3:
      pts( 0, 0) = -0.9352948105455122; pts( 0, 1) = -0.9352948105455122; pts( 0, 2) =  0.8058844316365366; 
      pts( 1, 0) =  0.8058844316365366; pts( 1, 1) = -0.9352948105455122; pts( 1, 2) = -0.9352948105455122; 
      pts( 2, 0) = -0.9352948105455122; pts( 2, 1) =  0.8058844316365366; pts( 2, 2) = -0.9352948105455122; 
      pts( 3, 0) = -0.9352948105455122; pts( 3, 1) = -0.9352948105455122; pts( 3, 2) = -0.9352948105455122; 
      pts( 4, 0) = -0.3804613914542760; pts( 4, 1) = -0.3804613914542760; pts( 4, 2) = -0.8586158256371720; 
      pts( 5, 0) = -0.8586158256371720; pts( 5, 1) = -0.3804613914542760; pts( 5, 2) = -0.3804613914542760; 
      pts( 6, 0) = -0.3804613914542760; pts( 6, 1) = -0.8586158256371720; pts( 6, 2) = -0.3804613914542760; 
      pts( 7, 0) = -0.3804613914542760; pts( 7, 1) = -0.3804613914542760; pts( 7, 2) = -0.3804613914542760; 
      pts( 8, 0) =  0.2331930661238736; pts( 8, 1) = -0.8792791169497158; pts( 8, 2) = -0.4746348322244420; 
      pts( 9, 0) =  0.2331930661238736; pts( 9, 1) = -0.8792791169497158; pts( 9, 2) = -0.8792791169497158; 
      pts(10, 0) = -0.8792791169497158; pts(10, 1) = -0.8792791169497158; pts(10, 2) =  0.2331930661238736; 
      pts(11, 0) = -0.4746348322244420; pts(11, 1) =  0.2331930661238736; pts(11, 2) = -0.8792791169497158; 
      pts(12, 0) = -0.8792791169497158; pts(12, 1) = -0.4746348322244420; pts(12, 2) =  0.2331930661238736; 
      pts(13, 0) = -0.8792791169497158; pts(13, 1) =  0.2331930661238736; pts(13, 2) = -0.8792791169497158; 
      pts(14, 0) = -0.4746348322244420; pts(14, 1) = -0.8792791169497158; pts(14, 2) =  0.2331930661238736; 
      pts(15, 0) = -0.8792791169497158; pts(15, 1) = -0.4746348322244420; pts(15, 2) = -0.8792791169497158; 
      pts(16, 0) = -0.8792791169497158; pts(16, 1) = -0.8792791169497158; pts(16, 2) = -0.4746348322244420; 
      pts(17, 0) = -0.8792791169497158; pts(17, 1) =  0.2331930661238736; pts(17, 2) = -0.4746348322244420; 
      pts(18, 0) = -0.4746348322244420; pts(18, 1) = -0.8792791169497158; pts(18, 2) = -0.8792791169497158; 
      pts(19, 0) =  0.2331930661238736; pts(19, 1) = -0.4746348322244420; pts(19, 2) = -0.8792791169497158; 
      break;

    case 4:
      pts( 0, 0) = -0.5000000000000000; pts( 0, 1) = -0.5000000000000000; pts( 0, 2) = -0.5000000000000000; 
      pts( 1, 0) = -0.9465264488912530; pts( 1, 1) = -0.9465264488912530; pts( 1, 2) =  0.8395793466737590; 
      pts( 2, 0) =  0.8395793466737590; pts( 2, 1) = -0.9465264488912530; pts( 2, 2) = -0.9465264488912530; 
      pts( 3, 0) = -0.9465264488912530; pts( 3, 1) =  0.8395793466737590; pts( 3, 2) = -0.9465264488912530; 
      pts( 4, 0) = -0.9465264488912530; pts( 4, 1) = -0.9465264488912530; pts( 4, 2) = -0.9465264488912530; 
      pts( 5, 0) = -0.9095091999689656; pts( 5, 1) = -0.0904908000310344; pts( 5, 2) = -0.9095091999689656; 
      pts( 6, 0) = -0.0904908000310344; pts( 6, 1) = -0.9095091999689656; pts( 6, 2) = -0.9095091999689656; 
      pts( 7, 0) = -0.9095091999689656; pts( 7, 1) = -0.0904908000310344; pts( 7, 2) = -0.0904908000310344; 
      pts( 8, 0) = -0.0904908000310344; pts( 8, 1) = -0.9095091999689656; pts( 8, 2) = -0.0904908000310344; 
      pts( 9, 0) = -0.0904908000310344; pts( 9, 1) = -0.0904908000310344; pts( 9, 2) = -0.9095091999689656; 
      pts(10, 0) = -0.9095091999689656; pts(10, 1) = -0.9095091999689656; pts(10, 2) = -0.0904908000310344; 
      pts(11, 0) =  0.4955197769636168; pts(11, 1) = -0.9217955187287024; pts(11, 2) = -0.6519287395062120; 
      pts(12, 0) =  0.4955197769636168; pts(12, 1) = -0.9217955187287024; pts(12, 2) = -0.9217955187287024; 
      pts(13, 0) = -0.9217955187287024; pts(13, 1) = -0.9217955187287024; pts(13, 2) =  0.4955197769636168; 
      pts(14, 0) = -0.6519287395062120; pts(14, 1) =  0.4955197769636168; pts(14, 2) = -0.9217955187287024; 
      pts(15, 0) = -0.9217955187287024; pts(15, 1) = -0.6519287395062120; pts(15, 2) =  0.4955197769636168; 
      pts(16, 0) = -0.9217955187287024; pts(16, 1) =  0.4955197769636168; pts(16, 2) = -0.9217955187287024; 
      pts(17, 0) = -0.6519287395062120; pts(17, 1) = -0.9217955187287024; pts(17, 2) =  0.4955197769636168; 
      pts(18, 0) = -0.9217955187287024; pts(18, 1) = -0.6519287395062120; pts(18, 2) = -0.9217955187287024; 
      pts(19, 0) = -0.9217955187287024; pts(19, 1) = -0.9217955187287024; pts(19, 2) = -0.6519287395062120; 
      pts(20, 0) = -0.9217955187287024; pts(20, 1) =  0.4955197769636168; pts(20, 2) = -0.6519287395062120; 
      pts(21, 0) = -0.6519287395062120; pts(21, 1) = -0.9217955187287024; pts(21, 2) = -0.9217955187287024; 
      pts(22, 0) =  0.4955197769636168; pts(22, 1) = -0.6519287395062120; pts(22, 2) = -0.9217955187287024; 
      pts(23, 0) =  0.0062372900291960; pts(23, 1) = -0.5535979240753700; pts(23, 2) = -0.8990414418784560; 
      pts(24, 0) =  0.0062372900291960; pts(24, 1) = -0.5535979240753700; pts(24, 2) = -0.5535979240753700; 
      pts(25, 0) = -0.5535979240753700; pts(25, 1) = -0.5535979240753700; pts(25, 2) =  0.0062372900291960; 
      pts(26, 0) = -0.8990414418784560; pts(26, 1) =  0.0062372900291960; pts(26, 2) = -0.5535979240753700; 
      pts(27, 0) = -0.5535979240753700; pts(27, 1) = -0.8990414418784560; pts(27, 2) =  0.0062372900291960; 
      pts(28, 0) = -0.5535979240753700; pts(28, 1) =  0.0062372900291960; pts(28, 2) = -0.5535979240753700; 
      pts(29, 0) = -0.8990414418784560; pts(29, 1) = -0.5535979240753700; pts(29, 2) =  0.0062372900291960; 
      pts(30, 0) = -0.5535979240753700; pts(30, 1) = -0.8990414418784560; pts(30, 2) = -0.5535979240753700; 
      pts(31, 0) = -0.5535979240753700; pts(31, 1) = -0.5535979240753700; pts(31, 2) = -0.8990414418784560; 
      pts(32, 0) = -0.5535979240753700; pts(32, 1) =  0.0062372900291960; pts(32, 2) = -0.8990414418784560; 
      pts(33, 0) = -0.8990414418784560; pts(33, 1) = -0.5535979240753700; pts(33, 2) = -0.5535979240753700; 
      pts(34, 0) =  0.0062372900291960; pts(34, 1) = -0.8990414418784560; pts(34, 2) = -0.5535979240753700; 
      break;

    case 5:
      pts( 0, 0) = -0.9700958696938816; pts( 0, 1) = -0.9700958696938816; pts( 0, 2) =  0.9102876090816449; 
      pts( 1, 0) =  0.9102876090816449; pts( 1, 1) = -0.9700958696938816; pts( 1, 2) = -0.9700958696938816; 
      pts( 2, 0) = -0.9700958696938816; pts( 2, 1) =  0.9102876090816449; pts( 2, 2) = -0.9700958696938816; 
      pts( 3, 0) = -0.9700958696938816; pts( 3, 1) = -0.9700958696938816; pts( 3, 2) = -0.9700958696938816; 
      pts( 4, 0) = -0.7310433304140120; pts( 4, 1) = -0.7310433304140120; pts( 4, 2) =  0.1931299912420360; 
      pts( 5, 0) =  0.1931299912420360; pts( 5, 1) = -0.7310433304140120; pts( 5, 2) = -0.7310433304140120; 
      pts( 6, 0) = -0.7310433304140120; pts( 6, 1) =  0.1931299912420360; pts( 6, 2) = -0.7310433304140120; 
      pts( 7, 0) = -0.7310433304140120; pts( 7, 1) = -0.7310433304140120; pts( 7, 2) = -0.7310433304140120; 
      pts( 8, 0) =  0.5599520168830801; pts( 8, 1) = -0.9318079576074770; pts( 8, 2) = -0.6963361016681260; 
      pts( 9, 0) =  0.5599520168830801; pts( 9, 1) = -0.9318079576074770; pts( 9, 2) = -0.9318079576074770; 
      pts(10, 0) = -0.9318079576074770; pts(10, 1) = -0.9318079576074770; pts(10, 2) =  0.5599520168830801; 
      pts(11, 0) = -0.6963361016681260; pts(11, 1) =  0.5599520168830801; pts(11, 2) = -0.9318079576074770; 
      pts(12, 0) = -0.9318079576074770; pts(12, 1) = -0.6963361016681260; pts(12, 2) =  0.5599520168830801; 
      pts(13, 0) = -0.9318079576074770; pts(13, 1) =  0.5599520168830801; pts(13, 2) = -0.9318079576074770; 
      pts(14, 0) = -0.6963361016681260; pts(14, 1) = -0.9318079576074770; pts(14, 2) =  0.5599520168830801; 
      pts(15, 0) = -0.9318079576074770; pts(15, 1) = -0.6963361016681260; pts(15, 2) = -0.9318079576074770; 
      pts(16, 0) = -0.9318079576074770; pts(16, 1) = -0.9318079576074770; pts(16, 2) = -0.6963361016681260; 
      pts(17, 0) = -0.9318079576074770; pts(17, 1) =  0.5599520168830801; pts(17, 2) = -0.6963361016681260; 
      pts(18, 0) = -0.6963361016681260; pts(18, 1) = -0.9318079576074770; pts(18, 2) = -0.9318079576074770; 
      pts(19, 0) =  0.5599520168830801; pts(19, 1) = -0.6963361016681260; pts(19, 2) = -0.9318079576074770; 
      pts(20, 0) =  0.1053112862120352; pts(20, 1) = -0.9075896991699967; pts(20, 2) = -0.2901318878720420; 
      pts(21, 0) =  0.1053112862120352; pts(21, 1) = -0.9075896991699967; pts(21, 2) = -0.9075896991699967; 
      pts(22, 0) = -0.9075896991699967; pts(22, 1) = -0.9075896991699967; pts(22, 2) =  0.1053112862120352; 
      pts(23, 0) = -0.2901318878720420; pts(23, 1) =  0.1053112862120352; pts(23, 2) = -0.9075896991699967; 
      pts(24, 0) = -0.9075896991699967; pts(24, 1) = -0.2901318878720420; pts(24, 2) =  0.1053112862120352; 
      pts(25, 0) = -0.9075896991699967; pts(25, 1) =  0.1053112862120352; pts(25, 2) = -0.9075896991699967; 
      pts(26, 0) = -0.2901318878720420; pts(26, 1) = -0.9075896991699967; pts(26, 2) =  0.1053112862120352; 
      pts(27, 0) = -0.9075896991699967; pts(27, 1) = -0.2901318878720420; pts(27, 2) = -0.9075896991699967; 
      pts(28, 0) = -0.9075896991699967; pts(28, 1) = -0.9075896991699967; pts(28, 2) = -0.2901318878720420; 
      pts(29, 0) = -0.9075896991699967; pts(29, 1) =  0.1053112862120352; pts(29, 2) = -0.2901318878720420; 
      pts(30, 0) = -0.2901318878720420; pts(30, 1) = -0.9075896991699967; pts(30, 2) = -0.9075896991699967; 
      pts(31, 0) =  0.1053112862120352; pts(31, 1) = -0.2901318878720420; pts(31, 2) = -0.9075896991699967; 
      pts(32, 0) =  0.0762086457760010; pts(32, 1) = -0.5436190778624780; pts(32, 2) = -0.9889704900510450; 
      pts(33, 0) =  0.0762086457760010; pts(33, 1) = -0.5436190778624780; pts(33, 2) = -0.5436190778624780; 
      pts(34, 0) = -0.5436190778624780; pts(34, 1) = -0.5436190778624780; pts(34, 2) =  0.0762086457760010; 
      pts(35, 0) = -0.9889704900510450; pts(35, 1) =  0.0762086457760010; pts(35, 2) = -0.5436190778624780; 
      pts(36, 0) = -0.5436190778624780; pts(36, 1) = -0.9889704900510450; pts(36, 2) =  0.0762086457760010; 
      pts(37, 0) = -0.5436190778624780; pts(37, 1) =  0.0762086457760010; pts(37, 2) = -0.5436190778624780; 
      pts(38, 0) = -0.9889704900510450; pts(38, 1) = -0.5436190778624780; pts(38, 2) =  0.0762086457760010; 
      pts(39, 0) = -0.5436190778624780; pts(39, 1) = -0.9889704900510450; pts(39, 2) = -0.5436190778624780; 
      pts(40, 0) = -0.5436190778624780; pts(40, 1) = -0.5436190778624780; pts(40, 2) = -0.9889704900510450; 
      pts(41, 0) = -0.5436190778624780; pts(41, 1) =  0.0762086457760010; pts(41, 2) = -0.9889704900510450; 
      pts(42, 0) = -0.9889704900510450; pts(42, 1) = -0.5436190778624780; pts(42, 2) = -0.5436190778624780; 
      pts(43, 0) =  0.0762086457760010; pts(43, 1) = -0.9889704900510450; pts(43, 2) = -0.5436190778624780; 
      pts(44, 0) = -0.6076324808508820; pts(44, 1) = -0.2953894798240120; pts(44, 2) = -0.8015885595010940; 
      pts(45, 0) = -0.6076324808508820; pts(45, 1) = -0.2953894798240120; pts(45, 2) = -0.2953894798240120; 
      pts(46, 0) = -0.2953894798240120; pts(46, 1) = -0.2953894798240120; pts(46, 2) = -0.6076324808508820; 
      pts(47, 0) = -0.8015885595010940; pts(47, 1) = -0.6076324808508820; pts(47, 2) = -0.2953894798240120; 
      pts(48, 0) = -0.2953894798240120; pts(48, 1) = -0.8015885595010940; pts(48, 2) = -0.6076324808508820; 
      pts(49, 0) = -0.2953894798240120; pts(49, 1) = -0.6076324808508820; pts(49, 2) = -0.2953894798240120; 
      pts(50, 0) = -0.8015885595010940; pts(50, 1) = -0.2953894798240120; pts(50, 2) = -0.6076324808508820; 
      pts(51, 0) = -0.2953894798240120; pts(51, 1) = -0.8015885595010940; pts(51, 2) = -0.2953894798240120; 
      pts(52, 0) = -0.2953894798240120; pts(52, 1) = -0.2953894798240120; pts(52, 2) = -0.8015885595010940; 
      pts(53, 0) = -0.2953894798240120; pts(53, 1) = -0.6076324808508820; pts(53, 2) = -0.8015885595010940; 
      pts(54, 0) = -0.8015885595010940; pts(54, 1) = -0.2953894798240120; pts(54, 2) = -0.2953894798240120; 
      pts(55, 0) = -0.6076324808508820; pts(55, 1) = -0.8015885595010940; pts(55, 2) = -0.2953894798240120; 
      break;

    case 6:
      pts( 0, 0) = -0.9462430511703660; pts( 0, 1) = -0.9462430511703660; pts( 0, 2) =  0.8387291535110980; 
      pts( 1, 0) =  0.8387291535110980; pts( 1, 1) = -0.9462430511703660; pts( 1, 2) = -0.9462430511703660; 
      pts( 2, 0) = -0.9462430511703660; pts( 2, 1) =  0.8387291535110980; pts( 2, 2) = -0.9462430511703660; 
      pts( 3, 0) = -0.9462430511703660; pts( 3, 1) = -0.9462430511703660; pts( 3, 2) = -0.9462430511703660; 
      pts( 4, 0) = -0.6257186483930600; pts( 4, 1) = -0.6257186483930600; pts( 4, 2) = -0.1228440548208200; 
      pts( 5, 0) = -0.1228440548208200; pts( 5, 1) = -0.6257186483930600; pts( 5, 2) = -0.6257186483930600; 
      pts( 6, 0) = -0.6257186483930600; pts( 6, 1) = -0.1228440548208200; pts( 6, 2) = -0.6257186483930600; 
      pts( 7, 0) = -0.6257186483930600; pts( 7, 1) = -0.6257186483930600; pts( 7, 2) = -0.6257186483930600; 
      pts( 8, 0) = -0.3557771363382860; pts( 8, 1) = -0.3557771363382860; pts( 8, 2) = -0.9326685909851420; 
      pts( 9, 0) = -0.9326685909851420; pts( 9, 1) = -0.3557771363382860; pts( 9, 2) = -0.3557771363382860; 
      pts(10, 0) = -0.3557771363382860; pts(10, 1) = -0.9326685909851420; pts(10, 2) = -0.3557771363382860; 
      pts(11, 0) = -0.3557771363382860; pts(11, 1) = -0.3557771363382860; pts(11, 2) = -0.3557771363382860; 
      pts(12, 0) = -0.9471516702558740; pts(12, 1) = -0.0528483297441260; pts(12, 2) = -0.9471516702558740; 
      pts(13, 0) = -0.0528483297441260; pts(13, 1) = -0.9471516702558740; pts(13, 2) = -0.9471516702558740; 
      pts(14, 0) = -0.9471516702558740; pts(14, 1) = -0.0528483297441260; pts(14, 2) = -0.0528483297441260; 
      pts(15, 0) = -0.0528483297441260; pts(15, 1) = -0.9471516702558740; pts(15, 2) = -0.0528483297441260; 
      pts(16, 0) = -0.0528483297441260; pts(16, 1) = -0.0528483297441260; pts(16, 2) = -0.9471516702558740; 
      pts(17, 0) = -0.9471516702558740; pts(17, 1) = -0.9471516702558740; pts(17, 2) = -0.0528483297441260; 
      pts(18, 0) = -0.7040905240547120; pts(18, 1) = -0.2959094759452880; pts(18, 2) = -0.7040905240547120; 
      pts(19, 0) = -0.2959094759452880; pts(19, 1) = -0.7040905240547120; pts(19, 2) = -0.7040905240547120; 
      pts(20, 0) = -0.7040905240547120; pts(20, 1) = -0.2959094759452880; pts(20, 2) = -0.2959094759452880; 
      pts(21, 0) = -0.2959094759452880; pts(21, 1) = -0.7040905240547120; pts(21, 2) = -0.2959094759452880; 
      pts(22, 0) = -0.2959094759452880; pts(22, 1) = -0.2959094759452880; pts(22, 2) = -0.7040905240547120; 
      pts(23, 0) = -0.7040905240547120; pts(23, 1) = -0.7040905240547120; pts(23, 2) = -0.2959094759452880; 
      pts(24, 0) =  0.4646198193858960; pts(24, 1) = -0.9580931155598880; pts(24, 2) = -0.5484335882661200; 
      pts(25, 0) =  0.4646198193858960; pts(25, 1) = -0.9580931155598880; pts(25, 2) = -0.9580931155598880; 
      pts(26, 0) = -0.9580931155598880; pts(26, 1) = -0.9580931155598880; pts(26, 2) =  0.4646198193858960; 
      pts(27, 0) = -0.5484335882661200; pts(27, 1) =  0.4646198193858960; pts(27, 2) = -0.9580931155598880; 
      pts(28, 0) = -0.9580931155598880; pts(28, 1) = -0.5484335882661200; pts(28, 2) =  0.4646198193858960; 
      pts(29, 0) = -0.9580931155598880; pts(29, 1) =  0.4646198193858960; pts(29, 2) = -0.9580931155598880; 
      pts(30, 0) = -0.5484335882661200; pts(30, 1) = -0.9580931155598880; pts(30, 2) =  0.4646198193858960; 
      pts(31, 0) = -0.9580931155598880; pts(31, 1) = -0.5484335882661200; pts(31, 2) = -0.9580931155598880; 
      pts(32, 0) = -0.9580931155598880; pts(32, 1) = -0.9580931155598880; pts(32, 2) = -0.5484335882661200; 
      pts(33, 0) = -0.9580931155598880; pts(33, 1) =  0.4646198193858960; pts(33, 2) = -0.5484335882661200; 
      pts(34, 0) = -0.5484335882661200; pts(34, 1) = -0.9580931155598880; pts(34, 2) = -0.9580931155598880; 
      pts(35, 0) =  0.4646198193858960; pts(35, 1) = -0.5484335882661200; pts(35, 2) = -0.9580931155598880; 
      pts(36, 0) =  0.2951151881739520; pts(36, 1) = -0.8060205337530680; pts(36, 2) = -0.6830741206678160; 
      pts(37, 0) =  0.2951151881739520; pts(37, 1) = -0.8060205337530680; pts(37, 2) = -0.8060205337530680; 
      pts(38, 0) = -0.8060205337530680; pts(38, 1) = -0.8060205337530680; pts(38, 2) =  0.2951151881739520; 
      pts(39, 0) = -0.6830741206678160; pts(39, 1) =  0.2951151881739520; pts(39, 2) = -0.8060205337530680; 
      pts(40, 0) = -0.8060205337530680; pts(40, 1) = -0.6830741206678160; pts(40, 2) =  0.2951151881739520; 
      pts(41, 0) = -0.8060205337530680; pts(41, 1) =  0.2951151881739520; pts(41, 2) = -0.8060205337530680; 
      pts(42, 0) = -0.6830741206678160; pts(42, 1) = -0.8060205337530680; pts(42, 2) =  0.2951151881739520; 
      pts(43, 0) = -0.8060205337530680; pts(43, 1) = -0.6830741206678160; pts(43, 2) = -0.8060205337530680; 
      pts(44, 0) = -0.8060205337530680; pts(44, 1) = -0.8060205337530680; pts(44, 2) = -0.6830741206678160; 
      pts(45, 0) = -0.8060205337530680; pts(45, 1) =  0.2951151881739520; pts(45, 2) = -0.6830741206678160; 
      pts(46, 0) = -0.6830741206678160; pts(46, 1) = -0.8060205337530680; pts(46, 2) = -0.8060205337530680; 
      pts(47, 0) =  0.2951151881739520; pts(47, 1) = -0.6830741206678160; pts(47, 2) = -0.8060205337530680; 
      pts(48, 0) =  0.5858785129392360; pts(48, 1) = -0.8047836742191160; pts(48, 2) = -0.9763111645010040; 
      pts(49, 0) =  0.5858785129392360; pts(49, 1) = -0.8047836742191160; pts(49, 2) = -0.8047836742191160; 
      pts(50, 0) = -0.8047836742191160; pts(50, 1) = -0.8047836742191160; pts(50, 2) =  0.5858785129392360; 
      pts(51, 0) = -0.9763111645010040; pts(51, 1) =  0.5858785129392360; pts(51, 2) = -0.8047836742191160; 
      pts(52, 0) = -0.8047836742191160; pts(52, 1) = -0.9763111645010040; pts(52, 2) =  0.5858785129392360; 
      pts(53, 0) = -0.8047836742191160; pts(53, 1) =  0.5858785129392360; pts(53, 2) = -0.8047836742191160; 
      pts(54, 0) = -0.9763111645010040; pts(54, 1) = -0.8047836742191160; pts(54, 2) =  0.5858785129392360; 
      pts(55, 0) = -0.8047836742191160; pts(55, 1) = -0.9763111645010040; pts(55, 2) = -0.8047836742191160; 
      pts(56, 0) = -0.8047836742191160; pts(56, 1) = -0.8047836742191160; pts(56, 2) = -0.9763111645010040; 
      pts(57, 0) = -0.8047836742191160; pts(57, 1) =  0.5858785129392360; pts(57, 2) = -0.9763111645010040; 
      pts(58, 0) = -0.9763111645010040; pts(58, 1) = -0.8047836742191160; pts(58, 2) = -0.8047836742191160; 
      pts(59, 0) =  0.5858785129392360; pts(59, 1) = -0.9763111645010040; pts(59, 2) = -0.8047836742191160; 
      pts(60, 0) =  0.0823688256004740; pts(60, 1) = -0.7328836785928640; pts(60, 2) = -0.4069979589137520; 
      pts(61, 0) = -0.7328836785928640; pts(61, 1) = -0.4069979589137520; pts(61, 2) =  0.0823688256004740; 
      pts(62, 0) = -0.7328836785928640; pts(62, 1) =  0.0823688256004740; pts(62, 2) = -0.4069979589137520; 
      pts(63, 0) = -0.4069979589137520; pts(63, 1) =  0.0823688256004740; pts(63, 2) = -0.9424871880938580; 
      pts(64, 0) = -0.4069979589137520; pts(64, 1) =  0.0823688256004740; pts(64, 2) = -0.7328836785928640; 
      pts(65, 0) = -0.9424871880938580; pts(65, 1) =  0.0823688256004740; pts(65, 2) = -0.4069979589137520; 
      pts(66, 0) = -0.4069979589137520; pts(66, 1) = -0.7328836785928640; pts(66, 2) = -0.9424871880938580; 
      pts(67, 0) = -0.4069979589137520; pts(67, 1) = -0.9424871880938580; pts(67, 2) =  0.0823688256004740; 
      pts(68, 0) = -0.9424871880938580; pts(68, 1) = -0.7328836785928640; pts(68, 2) = -0.4069979589137520; 
      pts(69, 0) =  0.0823688256004740; pts(69, 1) = -0.7328836785928640; pts(69, 2) = -0.9424871880938580; 
      pts(70, 0) = -0.9424871880938580; pts(70, 1) =  0.0823688256004740; pts(70, 2) = -0.7328836785928640; 
      pts(71, 0) = -0.4069979589137520; pts(71, 1) = -0.7328836785928640; pts(71, 2) =  0.0823688256004740; 
      pts(72, 0) = -0.9424871880938580; pts(72, 1) = -0.7328836785928640; pts(72, 2) =  0.0823688256004740; 
      pts(73, 0) =  0.0823688256004740; pts(73, 1) = -0.9424871880938580; pts(73, 2) = -0.4069979589137520; 
      pts(74, 0) = -0.4069979589137520; pts(74, 1) = -0.9424871880938580; pts(74, 2) = -0.7328836785928640; 
      pts(75, 0) = -0.7328836785928640; pts(75, 1) = -0.9424871880938580; pts(75, 2) = -0.4069979589137520; 
      pts(76, 0) =  0.0823688256004740; pts(76, 1) = -0.9424871880938580; pts(76, 2) = -0.7328836785928640; 
      pts(77, 0) = -0.7328836785928640; pts(77, 1) = -0.9424871880938580; pts(77, 2) =  0.0823688256004740; 
      pts(78, 0) = -0.7328836785928640; pts(78, 1) =  0.0823688256004740; pts(78, 2) = -0.9424871880938580; 
      pts(79, 0) = -0.9424871880938580; pts(79, 1) = -0.4069979589137520; pts(79, 2) =  0.0823688256004740; 
      pts(80, 0) =  0.0823688256004740; pts(80, 1) = -0.4069979589137520; pts(80, 2) = -0.9424871880938580; 
      pts(81, 0) =  0.0823688256004740; pts(81, 1) = -0.4069979589137520; pts(81, 2) = -0.7328836785928640; 
      pts(82, 0) = -0.7328836785928640; pts(82, 1) = -0.4069979589137520; pts(82, 2) = -0.9424871880938580; 
      pts(83, 0) = -0.9424871880938580; pts(83, 1) = -0.4069979589137520; pts(83, 2) = -0.7328836785928640; 
      break;

    default:
      ThrowException("WS_Tet_pts undefined for provided order!");
  }

  return pts;
}

mdvector<double> WS_Tet_weights(unsigned int P)
{
  unsigned int nPts = (P+1)*(P+2)*(P+3)/6;
  mdvector<double> weights({nPts});

  switch(P)
  {
    case 0:
      weights(0) = 2.0;
      break;

    case 1:
      weights( 0) =  0.3333333333333333;
      weights( 1) =  0.3333333333333333;
      weights( 2) =  0.3333333333333333;
      weights( 3) =  0.3333333333333333;
      break;

    case 2:
      weights( 0) =  0.0635108464576119;
      weights( 1) =  0.0635108464576119;
      weights( 2) =  0.0635108464576119;
      weights( 3) =  0.0635108464576119;
      weights( 4) =  0.1798816579171480;
      weights( 5) =  0.1798816579171480;
      weights( 6) =  0.1798816579171480;
      weights( 7) =  0.1798816579171480;
      weights( 8) =  0.1798816579171480;
      weights( 9) =  0.1798816579171480;
      break;

    case 3:
      weights( 0) =  0.0094227663926260;
      weights( 1) =  0.0094227663926260;
      weights( 2) =  0.0094227663926260;
      weights( 3) =  0.0094227663926260;
      weights( 4) =  0.1359158910531573;
      weights( 5) =  0.1359158910531573;
      weights( 6) =  0.1359158910531573;
      weights( 7) =  0.1359158910531573;
      weights( 8) =  0.0626648919625169;
      weights( 9) =  0.0626648919625169;
      weights(10) =  0.0626648919625169;
      weights(11) =  0.0626648919625169;
      weights(12) =  0.0626648919625169;
      weights(13) =  0.0626648919625169;
      weights(14) =  0.0626648919625169;
      weights(15) =  0.0626648919625169;
      weights(16) =  0.0626648919625169;
      weights(17) =  0.0626648919625169;
      weights(18) =  0.0626648919625169;
      weights(19) =  0.0626648919625169;
      break;

    case 4:
      weights( 0) =  0.1242327641593787;
      weights( 1) =  0.0029200618620517;
      weights( 2) =  0.0029200618620517;
      weights( 3) =  0.0029200618620517;
      weights( 4) =  0.0029200618620517;
      weights( 5) =  0.0333740527582328;
      weights( 6) =  0.0333740527582328;
      weights( 7) =  0.0333740527582328;
      weights( 8) =  0.0333740527582328;
      weights( 9) =  0.0333740527582328;
      weights(10) =  0.0333740527582328;
      weights(11) =  0.0191194226903553;
      weights(12) =  0.0191194226903553;
      weights(13) =  0.0191194226903553;
      weights(14) =  0.0191194226903553;
      weights(15) =  0.0191194226903553;
      weights(16) =  0.0191194226903553;
      weights(17) =  0.0191194226903553;
      weights(18) =  0.0191194226903553;
      weights(19) =  0.0191194226903553;
      weights(20) =  0.0191194226903553;
      weights(21) =  0.0191194226903553;
      weights(22) =  0.0191194226903553;
      weights(23) =  0.0639785777410072;
      weights(24) =  0.0639785777410072;
      weights(25) =  0.0639785777410072;
      weights(26) =  0.0639785777410072;
      weights(27) =  0.0639785777410072;
      weights(28) =  0.0639785777410072;
      weights(29) =  0.0639785777410072;
      weights(30) =  0.0639785777410072;
      weights(31) =  0.0639785777410072;
      weights(32) =  0.0639785777410072;
      weights(33) =  0.0639785777410072;
      weights(34) =  0.0639785777410072;
      break;

    case 5:
      weights( 0) =  0.0013830816448187;
      weights( 1) =  0.0013830816448187;
      weights( 2) =  0.0013830816448187;
      weights( 3) =  0.0013830816448187;
      weights( 4) =  0.0488388488540144;
      weights( 5) =  0.0488388488540144;
      weights( 6) =  0.0488388488540144;
      weights( 7) =  0.0488388488540144;
      weights( 8) =  0.0128022193865973;
      weights( 9) =  0.0128022193865973;
      weights(10) =  0.0128022193865973;
      weights(11) =  0.0128022193865973;
      weights(12) =  0.0128022193865973;
      weights(13) =  0.0128022193865973;
      weights(14) =  0.0128022193865973;
      weights(15) =  0.0128022193865973;
      weights(16) =  0.0128022193865973;
      weights(17) =  0.0128022193865973;
      weights(18) =  0.0128022193865973;
      weights(19) =  0.0128022193865973;
      weights(20) =  0.0219325302397643;
      weights(21) =  0.0219325302397643;
      weights(22) =  0.0219325302397643;
      weights(23) =  0.0219325302397643;
      weights(24) =  0.0219325302397643;
      weights(25) =  0.0219325302397643;
      weights(26) =  0.0219325302397643;
      weights(27) =  0.0219325302397643;
      weights(28) =  0.0219325302397643;
      weights(29) =  0.0219325302397643;
      weights(30) =  0.0219325302397643;
      weights(31) =  0.0219325302397643;
      weights(32) =  0.0204997022017747;
      weights(33) =  0.0204997022017747;
      weights(34) =  0.0204997022017747;
      weights(35) =  0.0204997022017747;
      weights(36) =  0.0204997022017747;
      weights(37) =  0.0204997022017747;
      weights(38) =  0.0204997022017747;
      weights(39) =  0.0204997022017747;
      weights(40) =  0.0204997022017747;
      weights(41) =  0.0204997022017747;
      weights(42) =  0.0204997022017747;
      weights(43) =  0.0204997022017747;
      weights(44) =  0.0391360157833640;
      weights(45) =  0.0391360157833640;
      weights(46) =  0.0391360157833640;
      weights(47) =  0.0391360157833640;
      weights(48) =  0.0391360157833640;
      weights(49) =  0.0391360157833640;
      weights(50) =  0.0391360157833640;
      weights(51) =  0.0391360157833640;
      weights(52) =  0.0391360157833640;
      weights(53) =  0.0391360157833640;
      weights(54) =  0.0391360157833640;
      weights(55) =  0.0391360157833640;
      break;

    case 6:
      weights( 0) =  0.0028599135257547;
      weights( 1) =  0.0028599135257547;
      weights( 2) =  0.0028599135257547;
      weights( 3) =  0.0028599135257547;
      weights( 4) =  0.0277688555876920;
      weights( 5) =  0.0277688555876920;
      weights( 6) =  0.0277688555876920;
      weights( 7) =  0.0277688555876920;
      weights( 8) =  0.0306675755590480;
      weights( 9) =  0.0306675755590480;
      weights(10) =  0.0306675755590480;
      weights(11) =  0.0306675755590480;
      weights(12) =  0.0096135147526067;
      weights(13) =  0.0096135147526067;
      weights(14) =  0.0096135147526067;
      weights(15) =  0.0096135147526067;
      weights(16) =  0.0096135147526067;
      weights(17) =  0.0096135147526067;
      weights(18) =  0.0410652255462827;
      weights(19) =  0.0410652255462827;
      weights(20) =  0.0410652255462827;
      weights(21) =  0.0410652255462827;
      weights(22) =  0.0410652255462827;
      weights(23) =  0.0410652255462827;
      weights(24) =  0.0058104597518187;
      weights(25) =  0.0058104597518187;
      weights(26) =  0.0058104597518187;
      weights(27) =  0.0058104597518187;
      weights(28) =  0.0058104597518187;
      weights(29) =  0.0058104597518187;
      weights(30) =  0.0058104597518187;
      weights(31) =  0.0058104597518187;
      weights(32) =  0.0058104597518187;
      weights(33) =  0.0058104597518187;
      weights(34) =  0.0058104597518187;
      weights(35) =  0.0058104597518187;
      weights(36) =  0.0114580409037773;
      weights(37) =  0.0114580409037773;
      weights(38) =  0.0114580409037773;
      weights(39) =  0.0114580409037773;
      weights(40) =  0.0114580409037773;
      weights(41) =  0.0114580409037773;
      weights(42) =  0.0114580409037773;
      weights(43) =  0.0114580409037773;
      weights(44) =  0.0114580409037773;
      weights(45) =  0.0114580409037773;
      weights(46) =  0.0114580409037773;
      weights(47) =  0.0114580409037773;
      weights(48) =  0.0064840852065493;
      weights(49) =  0.0064840852065493;
      weights(50) =  0.0064840852065493;
      weights(51) =  0.0064840852065493;
      weights(52) =  0.0064840852065493;
      weights(53) =  0.0064840852065493;
      weights(54) =  0.0064840852065493;
      weights(55) =  0.0064840852065493;
      weights(56) =  0.0064840852065493;
      weights(57) =  0.0064840852065493;
      weights(58) =  0.0064840852065493;
      weights(59) =  0.0064840852065493;
      weights(60) =  0.0207935201043453;
      weights(61) =  0.0207935201043453;
      weights(62) =  0.0207935201043453;
      weights(63) =  0.0207935201043453;
      weights(64) =  0.0207935201043453;
      weights(65) =  0.0207935201043453;
      weights(66) =  0.0207935201043453;
      weights(67) =  0.0207935201043453;
      weights(68) =  0.0207935201043453;
      weights(69) =  0.0207935201043453;
      weights(70) =  0.0207935201043453;
      weights(71) =  0.0207935201043453;
      weights(72) =  0.0207935201043453;
      weights(73) =  0.0207935201043453;
      weights(74) =  0.0207935201043453;
      weights(75) =  0.0207935201043453;
      weights(76) =  0.0207935201043453;
      weights(77) =  0.0207935201043453;
      weights(78) =  0.0207935201043453;
      weights(79) =  0.0207935201043453;
      weights(80) =  0.0207935201043453;
      weights(81) =  0.0207935201043453;
      weights(82) =  0.0207935201043453;
      weights(83) =  0.0207935201043453;
      break;

    default:
      ThrowException("WS_Tet_weights undefined for provided order!");
  }

  return weights;
}
